Section: New Results

Inverse Problems

Reconstruction of an elastic scatterer immersed in a homogeneous fluid

Participants : Hélène Barucq, Rabia Djellouli, Élodie Estecahandy.

The determination of the shape of an obstacle from its effects on known acoustic or electromagnetic waves is an important problem in many technologies such as sonar, radar, geophysical exploration, medical imaging and nondestructive testing. This inverse obstacle problem (IOP) is difficult to solve, especially from a numerical viewpoint, because it is ill-posed and nonlinear. Its investigation requires as a prerequisite the fundamental understanding of the theory for the associated direct scattering problem, and the mastery of the corresponding numerical solution methods.

In this work, we are interested in retrieving the shape of an elastic obstacle from the knowledge of some scattered far-field patterns, and assuming certain characteristics of the surface of the obstacle. The corresponding direct elasto-acoustic scattering problem consists in the scattering of time-harmonic acoustic waves by an elastic obstacle ${\Omega }^s$ embedded in a homogeneous medium ${\Omega }^f$, that can be formulated as follows:

where $p$ is the fluid pressure in ${\Omega }^f$ whereas $u$ is the displacement field in ${\Omega }^s$, and $\sigma \left(u\right)$ represents the stress tensor of the elastic material.

This boundary value problem has been investigated mathematically and results pertaining to the existence, uniqueness and regularity can be found in [65] and the references therein, among others. We propose a solution methodology based on a regularized Newton-type method for solving the IOP. The proposed method is an extension of the regularized Newton algorithm developed for solving the case where only Helmholtz equation is involved, that is the acoustic case by impenetrable scatterers [55] . The direct elasto-acoustic scattering problem defines an operator $F:\Gamma \to p_{\infty }$ which maps the boundary $\Gamma$ of the scatterer ${\Omega }^s$ onto the far-field pattern $p_{\infty }$. Hence, given one or several measured far-field patterns ${\tilde{p}}_{\infty }\left(\widehat{x}\right)$, corresponding to one or several given directions $d$ and wavenumbers $k$, one can formulate IOPs as follows:

We propose a solution methodology based on a regularized Newton-type method to solve this inverse obstacle problem. At each Newton iteration, we solve the forward problem using a finite element solver based on discontinuous Galerkin approximations, and equipped with high-order absorbing boundary conditions. We have first characterized the Fréchet derivatives of the scattered field. They are solution to the same boundary value problem as the direct problem with other transmission conditions. This work has been presented both in FACM11 and in WAVES 2011. A paper has been submitted.

Seismic data interpretation using the Hough transform and principal component analysis

Participants : M-.G Orozco-del-Castillo, Carlos Ortiz-Aleman, Roland Martin, Rafael Avila-Carrera, Alejandro Rodriguez-Castellanos.

In [29] , two novel image processing techniques are applied to detect and delineate complex salt bodies from seismic exploration profiles: Hough transform and principal component analysis (PCA). It is well recognized by the geophysical community that the lack of resolution and poor structural identification in seismic data recorded at sub-salt plays represent severe technical and economical problems. Under such circumstances, seismic interpretation based only on the human-eye is inaccurate. Additionally, petroleum field development decisions and production planning depend on good-quality seismic images that generally are not feasible in salt tectonics areas. In spite of this, morphological erosion, region growing and, especially, a generalization of the Hough transform (closely related to the Radon transform) are applied to build parabolic shapes that are useful in the idealization and recognition of salt domes from 2D seismic profiles. In a similar way, PCA is also used to identify shapes associated with complex salt bodies in seismic profiles extracted from 3D seismic data. To show the validity of the new set of seismic results, comparisons between both image processing techniques are exhibited. It is remarkable that the main contribution of this work is oriented in providing the seismic interpreters with new semi-automatic computational tools. The novel image processing approaches presented here may be helpful in the identification of diapirs and other complex geological features from seismic images. Conceivably, in the near future, a new branch of seismic attributes could be recognized by geoscientists and engineers based on the encouraging results reported here.

Gravimetry Inversion

Participants : Roland Martin, Dimitri Komatitsch, Mark Jessel, Stéphane Perrouty, Vadim Monteiller.

In order to improve the subsoil images in regions which are not well covered by a dense seismic array or can not be well retrieved by using seismic imaging techniques alone (salty dome regions like in the Gulf of Mexico for instance), we have been developping new gravity imaging techniques using supercomputing. In regions like Ghana or the Chicxulub crater located in Yucatan plate (Mexico), 3D sensitivity kernels are calculated for gravity potential data sets measured over 2000 up to 10000 locations randomly distributed in space. The density anomaly computational domain covers a 250 km $\times$ 250 km $\times$ 20 km volume in these two regions. For instance for the Ghana region two resolutions are taken : 1 point each 2 kms in the horizontal plane and 200 m in the vertical direction in the less accurate configuration and one point each 500 m in the horizontal plane and one point each 100 m in the vertical direction for the most accurate configuration. The gravity anomalies are inverted using an optimized least-square method applied to a sensitivity kernel of ${10}^{1}0$ up to $4\times {10}^{1}1$ elements . The least-square method using a $L^2$-norm or $L^1$-norm has been implemented on hybrid multi-CPU/multi-GPU Titane machines at CCRT of the French Nuclear Energy Agency. $L^1$ norm gives us sharper boundaries of the density structures when compared to $L^2$-norm solutions but these $L^1$ solutions are obtained at the expense of one order of magnitude in the number of iterations necessary to the inversion. The $L^2$-norm gives slightly smoother solutions but is much more faster than $L^1$ norm by at least one order of magnitude in terms of acceleration. The optimized CPU version has allowed us to reduce drastically the computation time form 5 hours on 512 processors to 25 minutes. Furtermore the multi-GPU version has decreased this computational time around 15 minutes. Our collaboration with geologists of IRD (Mark Jessell and Stephane Perrouty) has allowed us to determine a reallistic a priori model of Ghana with different resolutions in order to obtain more reallistic models after inversion. We are still optimizing the multi-GPU code, with the challenging goal of obtaining results lower than 10 minutes and then obtaining an acceleration factor of at least 4 when compared to the optimized multi-CPU iunversion code. By now the code is further optimized and more improvements are alreday in curse in terms of better preconditionning, automatic multi-resolution procedure implementation, gravimetry-seismic joint inversion and extension to the global earth imaging. An international article in a scientific peer-reviewed journal is in preparation.